for (i in 1:length(params))
  print(paste('Parameter:', names(params)[i], ' - Value:', params[[i]], '- Class:', class(params[[i]])))
## [1] "Parameter: Dataset  - Value: CHD2_iPSCs_and_organoids_PublicRepo - Class: character"
## [1] "Parameter: CountFile  - Value: /group/testa/Project/CHD2/BulkRNAseq/data/Counts/CHD2_bulkRNAseqMergedCounts.txt - Class: character"
## [1] "Parameter: DesignFile  - Value: /group/testa/Project/CHD2/BulkRNAseq/data/SampleSheets/CHD2_BulkRNAseq_mergedRuns_SampleSheet.txt - Class: character"
## [1] "Parameter: GeneAnnotationFile  - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/0a.GeneAnnotation/Output/AnnotationEns101.txt - Class: character"
## [1] "Parameter: InputFolder  - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/1.Exploratory/Input/ - Class: character"
## [1] "Parameter: SavingFolder  - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/1.Exploratory/Output/Savings/ - Class: character"
## [1] "Parameter: OutputFolder  - Value: /group/testa/Project/CHD2/BulkRNAseq/results/PublicRepo/1.Exploratory/Output/Figures/ - Class: character"
## [1] "Parameter: CpmFilt  - Value: 2 - Class: integer"
## [1] "Parameter: SampleFilt  - Value: 3 - Class: integer"

Parameters specified while running this markdown:

Dataset: CHD2_iPSCs_and_organoids_PublicRepo
CountFile: countsmatrix.txt. Input files with read counts aggregated at the gene level.
DesignFile: CHD2_BulkRNAseq_mergedRuns_SampleSheet.txt
GeneAnnotationFile: reporting gene metadata. It must contain a ‘Gene’ column reporting GenecodeIDs and and ‘EnsGene’ column with Ensembl Gene ID.
CpmFilt: threshold in cpm (count per million) for the filtering of expressed genes. Set = 2.
SampleFilt: minimum number of samples that have to respect the threshold on cpm for the filtering of expressed genes. Set = 3

Libraries

library(DT)
library(gridExtra)
library(SummarizedExperiment)
library(ggplot2)
library(dplyr)
library(edgeR)
library(org.Hs.eg.db)
library(RNASeqBulkExploratory)
source("/group/testa/Project/CHD2/BulkRNAseq/data/HelperFunctions/HelperFunctions.R")
source("/group/testa/Project/CHD2/BulkRNAseq/data/HelperFunctions/PCA_helper_OL.R")
source("/group/testa/Project/CHD2/BulkRNAseq/data/HelperFunctions/loadingsPCAZeroVar.R")

Dataset <- params$Dataset
OutputFolder <- params$OutputFolder
SavingFolder <- params$SavingFolder

OutputFolder <- ifelse(is.null(params$OutputFolder), getwd(), params$OutputFolder) 
if (dir.exists(OutputFolder) == FALSE) {
  dir.create(OutputFolder, recursive=FALSE)
}

SavingFolder <- ifelse(is.null(params$SavingFolder), getwd(), params$SavingFolder) 
if (dir.exists(SavingFolder) == FALSE) {
  dir.create(SavingFolder, recursive=FALSE)
}

All samples: iPSCs, Day25 CO, Day50 CO

2. Data Upload

2.1 Read Count Matrix

RawCounts <- read.table(params$CountFile, sep='\t', header=T)

rownames(RawCounts) <- NULL
RawCounts <- RawCounts %>% dplyr::select(Gene, everything())

colSums(RawCounts[, -1])
##  OLE_001  OLE_002  OLE_004  OLE_005  OLE_007  OLE_008  OLE_010  OLE_011  OLE_013  OLE_014 
## 19060905 19766540 16841865 22320983 17458719 19168357 14691167 29931179 20824248 42056466 
##  OLE_016  OLE_017  OLE_026  OLE_027  OLE_028  OLE_029  OLE_030  OLE_031  OLE_032  OLE_033 
## 15707839 17807937 23426758 21638096 24224315 23238830 21819104 24407378 23196208 21580829 
##  OLE_034  OLE_035  OLE_036  OLE_037  OLE_038  OLE_039  OLE_040  OLE_041  OLE_042  OLE_043 
## 23510605 18957844 21991122 22781664 21287851 19458773 19274390 23971783 21355449 22917288 
##  OLE_044  OLE_045  OLE_046  OLE_047  OLE_048  OLE_049 
## 20795924 19222756 13134124 25570596 18118198 20488353

libsizeplot <- data.frame(libsize = colSums(RawCounts[, -1])) %>% 
  ggplot(., aes(y=rownames(.), x=libsize)) +
    geom_col(size=0.75, col='black') +
    labs(x='Library Size', y='Sample') +
    theme_bw() + 
    theme(axis.text = element_text(size=13, colour = 'black'), legend.key.size = unit(2.5, "line"), 
        legend.text = element_text(size=15), axis.title = element_text(size=15))
                       

libsizeplot

length(which(duplicated(RawCounts$Gene) == TRUE))
## [1] 0
length(RawCounts$Gene)
## [1] 60237
length(unique(RawCounts$Gene))
## [1] 60237

2.2 Retrieve Annotation from GeneAnnotation table and match Gene columns

GeneAnnotation <- read.table(params$GeneAnnotationFile, sep='\t', header=TRUE)

GeneAnnotation <- GeneAnnotation[match(RawCounts$Gene, GeneAnnotation$Gene), ] #give same order of dataset

if(!identical(nrow(RawCounts), nrow(GeneAnnotation))){
  stop('RawCounts and GeneAnnotation have different number of rows')
}

if(all(!identical(RawCounts$Gene, GeneAnnotation$Gene))){
  stop('Genes in RawCounts and GeneAnnotation are different')
}

dim(RawCounts)
## [1] 60237    37
dim(GeneAnnotation)
## [1] 60237     9

identical(RawCounts$Gene, GeneAnnotation$Gene)
## [1] TRUE

2.3 GeneAnnotation table sanity checks

if ('hgnc_symbol' %in% colnames(GeneAnnotation)){
  GeneAnnotation <- dplyr::rename(GeneAnnotation, HGNCSymbol=hgnc_symbol)
}

if ('external_gene_name' %in% colnames(GeneAnnotation)){
  GeneAnnotation <- dplyr::rename(GeneAnnotation, GeneName=external_gene_name)
}

if ('gene_biotype' %in% colnames(GeneAnnotation)){
  GeneAnnotation <- dplyr::rename(GeneAnnotation, GeneBiotype=gene_biotype)
}

if ('chromosome_name' %in% colnames(GeneAnnotation)){
  GeneAnnotation <- dplyr::rename(GeneAnnotation, Chr=chromosome_name)
}

if ('start_position' %in% colnames(GeneAnnotation)){
  GeneAnnotation <- dplyr::rename(GeneAnnotation, Start=start_position)
}

if ('end_position' %in% colnames(GeneAnnotation)){
  GeneAnnotation <- dplyr::rename(GeneAnnotation, End=end_position)
}

2.4 Specify experimental design and generate sample table

2.4.1 Import SampleSheet

Design <- read.table(params$DesignFile, sep='\t', header=TRUE)

2.4.2 Interactive table

as.data.frame(lapply(Design, factor)) %>%
  DT::datatable(class='hover', rownames=FALSE, caption='Sample identity and attributes',
                filter='top', escape=TRUE, extension='Buttons',
                options=list(pageLength=11, scrollX=T, dom='Bfrtip', buttons=list(I('colvis'), c('csv', 'excel'))))

2.5 Set ENSG as rownames and drop gene columns (but keep for later use)

if(all(RawCounts$Gene %in% GeneAnnotation$Gene) == FALSE){
    stop('ERROR! Count matrix rownames and annotation matrix are not consistent.') #final check that dataset and annotation are consistent
}

dim(RawCounts)
## [1] 60237    37
dim(GeneAnnotation)
## [1] 60237     9

3. Read Counts density

DensityRaw <- readCountDensityplot(countTable=RawCounts, GeneCol ='Gene', dgelist=FALSE, #Change Gene to gene to match counts matrix colname
                                   plotTitle='Raw Counts before filtering', prior.count=0.25, adjustment=0.5)
DensityRaw

4. SummarizedExperiment

4.1 Consistency checks

tmp <- RawCounts
o <- GeneAnnotation

#CountsMatrix
rownames(tmp) <- tmp$Gene
checkCounts(tmp, showDuplicates = FALSE)
## No duplicated genes to remove

#Design
checkDesign(Design, tmp[,-1]) #without the Gene column (in position 1)

#Annotation
o$EnsGene <- o$Gene #overwrite EnsGene column with info from Gene (containing gene version)
checkGeneAnno(o, tmp[,-1]) #without the Gene column (in position 1)

rm(tmp)
rm(o)

4.2 Create SummatizedExperiment object (direct function)

SE <- createSE(RawCounts, Design, GeneAnnotation, showDuplicates = TRUE, roundCounts = TRUE)
## No duplicated genes to remove

#colData(SE)

4.3 Set ENSG with version as rownames

rownames(SE) <- rowData(SE)$Gene #ENSG.version are unique

5. Biotypes selection

The purpose of this step is to discard all RNA species except for coding genes and lncRNAs.

5.1 Biotype-selected SE

SE_Bio <- biotypeSelectSE(SE, lncRNA=TRUE, otherBios=NULL, showRemoved=TRUE)
## UPDATED library.size slot
## A total of 23439 rows out of 60237 were discarded, belonging to the following biotypes:
## 
##                           artifact                          IG_C_gene 
##                                 19                                 14 
##                    IG_C_pseudogene                          IG_D_gene 
##                                  9                                 30 
##                          IG_J_gene                    IG_J_pseudogene 
##                                 18                                  3 
##                      IG_pseudogene                          IG_V_gene 
##                                  1                                144 
##                    IG_V_pseudogene                              miRNA 
##                                184                               1826 
##                           misc_RNA                            Mt_rRNA 
##                               2160                                  2 
##                            Mt_tRNA               processed_pseudogene 
##                                 22                              10115 
##                         pseudogene                           ribozyme 
##                                 15                                  8 
##                               rRNA                    rRNA_pseudogene 
##                                 26                                486 
##                             scaRNA                              scRNA 
##                                 48                                  1 
##                             snoRNA                              snRNA 
##                                924                               1836 
##                               sRNA                                TEC 
##                                  5                               1041 
##                          TR_C_gene                          TR_D_gene 
##                                  6                                  4 
##                          TR_J_gene                    TR_J_pseudogene 
##                                 78                                  4 
##                          TR_V_gene                    TR_V_pseudogene 
##                                106                                 33 
##   transcribed_processed_pseudogene     transcribed_unitary_pseudogene 
##                                508                                152 
## transcribed_unprocessed_pseudogene    translated_processed_pseudogene 
##                                953                                  2 
##  translated_unprocessed_pseudogene                 unitary_pseudogene 
##                                  3                                 98 
##             unprocessed_pseudogene                          vault_RNA 
##                               2554                                  1
## A total of 36607 rows out of 60237 were kept, belonging to the following biotypes:
## 
##         lncRNA protein_coding 
##          16735          19872
RemovedGenes <- SummarizedExperiment::rowData(SE)[!(SummarizedExperiment::rowData(SE) %in% 
                                                      SummarizedExperiment::rowData(SE_Bio)), ]

After the gene selection based on biotypes, the dataset is structured in 36607 genes measured in 36 samples. 23630 genes have been discarded.

5.2 Diagnostic Barplots

Percentage of Reads associated to most expressed gene

plotCompareMax(SE, SE_Bio, PlotTitle=NULL, Interactive=TRUE)
## Genes in the first dataset:  60237
## Genes in the second dataset:  36607

Percentage of Removed Reads after biotype selection

plotComparePerc(SE, SE_Bio, PlotTitle=NULL, Interactive=TRUE)
## Genes in the first dataset:  60237
## Genes in the second dataset:  36607

The dataset is structured in 60237 genes measured in 36 samples.

6. Gene expression filtering and SE normalization

6.1 Setting gene filtering thresholds

CpmFilt <- params$CpmFilt
if(is.numeric(CpmFilt)==FALSE){
  stop('ERROR! Cpm threshold for gene filtering has a non-numeric value')
}

SampleFilt <- ifelse(params$SampleFilt=='Default', dim(SE_Bio)[2]/2, params$SampleFilt)
if(is.numeric(SampleFilt)==FALSE){
  stop('ERROR! Sample threshold for gene filtering has a non-numeric value')
}

Genes having an expression lower than 2 CPM in at least 3 samples are filtered out.

6.2 Filter out low expressed genes

Compute temporary cpm (without TMM normalization, not stored in SE)
Filter out lowly-expressed genes

SE$lib.size
##  OLE_001  OLE_002  OLE_004  OLE_005  OLE_007  OLE_008  OLE_010  OLE_011  OLE_013  OLE_014 
## 19060754 19766345 16841681 22320810 17458566 19168180 14690968 29931004 20824033 42056306 
##  OLE_016  OLE_017  OLE_026  OLE_027  OLE_028  OLE_029  OLE_030  OLE_031  OLE_032  OLE_033 
## 15707707 17807771 23426589 21637962 24224175 23238642 21818959 24407222 23196008 21580680 
##  OLE_034  OLE_035  OLE_036  OLE_037  OLE_038  OLE_039  OLE_040  OLE_041  OLE_042  OLE_043 
## 23510421 18957702 21990976 22781464 21287710 19458674 19274220 23971621 21355321 22917142 
##  OLE_044  OLE_045  OLE_046  OLE_047  OLE_048  OLE_049 
## 20795731 19222606 13133999 25570386 18117992 20488146
SE_Filt <- filterSE(SE_Bio, cpmTh=CpmFilt, sampleTh=SampleFilt)
## UPDATED library.size slot
SE_Filt$lib.size
##  OLE_001  OLE_002  OLE_004  OLE_005  OLE_007  OLE_008  OLE_010  OLE_011  OLE_013  OLE_014 
## 18732742 19237611 16509621 21759751 17191349 18838759 14414057 29304725 20371687 41162813 
##  OLE_016  OLE_017  OLE_026  OLE_027  OLE_028  OLE_029  OLE_030  OLE_031  OLE_032  OLE_033 
## 15413332 17431247 23018916 21270865 23750876 22832395 21382838 24018734 22772845 21218735 
##  OLE_034  OLE_035  OLE_036  OLE_037  OLE_038  OLE_039  OLE_040  OLE_041  OLE_042  OLE_043 
## 23117391 18607665 21587839 22389686 20883465 19063501 18940436 23556871 20973750 22491864 
##  OLE_044  OLE_045  OLE_046  OLE_047  OLE_048  OLE_049 
## 20411058 18863938 12880577 25107868 17767668 20056391

6.3 Normalize

Re-calculate library size and calculate TMM
Calculate cpm and logCpm (with normalization) and store in SE_Filt

SE_Filt <- normTmmSE(SE_Filt, useNormFactors=TRUE, priorCount=0.25)  
## Warning in normTmmSE(SE_Filt, useNormFactors = TRUE, priorCount = 0.25): The following colData
## columns will be recomputed: lib.size

After filtering out lowly-expressed genes, the dataset is structured in 16225 genes measured in 36 samples.

7. Diagnostic barplots on normalized data

7.1 Barplot for library size

Barplot showing for each sample the initial library size; plot is centered on the mean value
Barplot with the same visualization after filtering on Cpm threshold

Library size before filtering

librarySizeBarplot(SE, PlotTitle='Before Filtering', Interactive=TRUE)

Library size after filtering

b1 <- librarySizeBarplot(SE, PlotTitle='Before Biotype Selection', Interactive=FALSE)
lim_y <- layer_scales(b1)$y$range$range

librarySizeBarplot(SE_Filt, PlotTitle='After Cpm Filtering', Interactive=TRUE, ylim=lim_y)

7.2 Barplot for TMM and expressed genes

Barplot showing for each sample the TMM value as calculated by edgeR and used in the normalization
Barplot represents in each sample the number of ‘expressed’ genes. The evaluation is done on SE before the gene filtering step and applying the specified Cpm threshold separately to each sample.

TMM Barplot

tMMBarplot(SE_Filt, Interactive=TRUE)

Barplot for genes with Cpm higher than threshold

expressedGeneBarplot(SE, PlotTitle='Before Filtering', cpmth=CpmFilt, Interactive=TRUE)

8. Density plot for Count distribution after filtering

Similarly to raw counts, I represent by density plot the count distribution also after filtering.

readCountDensityplot(assays(SE_Filt)$counts, plotTitle='Filtered Dataset')

9. Sample-to-sample correlation

Correlation matrix across samples calculated on the basis of the Spearman correlation. Heatmap annotation takes into consideration sample ‘Condition’.

# Set heatmap size on the basis of the number of samples
Width <- 7 + (dim(SE_Filt)[2]- 8) * 0.24
Width <- ifelse(Width<=14, Width, 14)
Height <- 5 + (dim(SE_Filt)[2]- 8) * 0.22

sampleCorrHeatmap(SE_Filt, cor_method='spearman', plotTitle='Filtered dataset', annotation_colors=NULL, display_numbers=TRUE, annotation_col = c('Timepoint', 'Genotype'))

## There are many samples. For readability "display_numbers" will be set to FALSE

10. Principal Component Analysis

if('Color' %in% names(Design)){
  Cols <- dplyr::select(Design, c(Genotype, Color)) %>% dplyr::arrange(Genotype) 
  Cols <- as.character(Cols[!duplicated(Cols$Genotype), ]$Color)
  } else {
    Cols <- NULL
    }

10.1 PCA with Features

All Samples

#PCA1 <- loadingsPCA(SE_Filt, components = c(1,2), condition = c('Genotype', 'Timepoint'), nFeatures = 10)

SE_Filt_GeneNames <- SE_Filt
SE_Filt_GeneNames <- SE_Filt_GeneNames[!duplicated(rowData(SE_Filt_GeneNames)$GeneName), ] #drop non-unique GeneNames
rownames(SE_Filt_GeneNames) <- rowData(SE_Filt_GeneNames)$GeneName #use GeneNames as rownames

loadingsPCA(SE_Filt_GeneNames, components = c(1,2), condition = c('Genotype', 'Timepoint'), nFeatures = 15)

Organoids

SE_Filt_GeneNames_organoids <- SE_Filt_GeneNames[, colData(SE_Filt_GeneNames)$Timepoint != 'd0']

loadingsPCAZeroVar(SE_Filt_GeneNames_organoids, components = c(1,2), condition = c('Genotype', 'Timepoint'), nFeatures = 15)

## Genes removed due to zero variance:
## [1] "DPEP3"     "ILDR1"     "CDA"       "GDF3"      "MFSD6L"    "TCL1B"     "LNCPRESS2"

ggsave(filename = paste0(OutputFolder, 'PCA_organoids_loadings_bulkRNAseq.pdf'), width = 5, height = 5)

Organoids day25

loadingsPCAZeroVar(SE_Filt_GeneNames_organoids[, colData(SE_Filt_GeneNames_organoids)$Timepoint == 'd25'], components = c(1,2), condition = c('Genotype', 'Timepoint'), nFeatures = 15)

## Genes removed due to zero variance:
##  [1] "SERPINB3"  "CDX4"      "MYBPH"     "DPEP3"     "ILDR1"     "UGT3A1"    "TRIM48"   
##  [8] "CDA"       "NEUROD6"   "CPLX4"     "KLK13"     "TRIM49B"   "GDF3"      "MFSD6L"   
## [15] "SERPINB4"  "TCL1B"     "LINC01108" "LNCPRESS2"

Organoids day50

loadingsPCAZeroVar(SE_Filt_GeneNames_organoids[, colData(SE_Filt_GeneNames_organoids)$Timepoint == 'd50'], components = c(1,2), condition = c('Genotype', 'Timepoint'), nFeatures = 15)

## Genes removed due to zero variance:
##  [1] "TMPRSS11E" "CR2"       "TNFSF11"   "OVOL2"     "RBBP8NL"   "DPEP3"     "ILDR1"    
##  [8] "CDA"       "DPPA2"     "TEX13B"    "TRIML2"    "GDF3"      "MFSD6L"    "LCN10"    
## [15] "CR1L"      "TCL1B"     "DANT1"     "LINC00992" "LNCPRESS2"

10.2 PCA minimal

pcaGenotype1 <- pcaSE(SE_Filt, PlotTitle='First-Second Component', condition="Genotype", Interactive=TRUE)

pcaGenotype2 <- pcaSE(SE_Filt, PlotTitle='Second-Third Component', components = c(2,3), condition="Genotype", Interactive=TRUE)

pcaTimepoint1 <- pcaSE(SE_Filt, PlotTitle='First-Second Component', condition="Timepoint", Interactive=TRUE)

pcaTimepoint2 <- pcaSE(SE_Filt, PlotTitle='Second-Third Component', components = c(2,3), condition="Timepoint", Interactive=TRUE)

pcaBatch <- pcaSE(SE_Filt, PlotTitle='First-Second Component', condition="Batch", Interactive=TRUE)

Genotype PC1-PC2

pcaGenotype1

Genotype PC2-PC3

pcaGenotype2

Timepoint PC1-PC2

pcaTimepoint1

Timepoint PC2-PC3

pcaTimepoint2

Batch PC1-PC2

pcaBatch

10.3 Timepoints-wise

#iPSCs
ips.se <- SE_Filt[, SE_Filt$Timepoint == 'd0']
pcaIPS <- pcaSE(ips.se, PlotTitle='First-Second Component iPSCs', condition="Genotype", Interactive=TRUE)

#Day25
d25.se <- SE_Filt[, SE_Filt$Timepoint == 'd25']
pcaD25 <- pcaSE(d25.se, PlotTitle='First-Second Component CO day25', condition="Genotype", Interactive=TRUE)

#Day50
d50.se <- SE_Filt[, SE_Filt$Timepoint == 'd50']
pcaD50 <- pcaSE(d50.se, PlotTitle='First-Second Component CO day50', condition="Genotype", Interactive=TRUE)

iPSCs

pcaIPS

Day25

pcaD25

Day50

pcaD50

10.3 Organoids

ORG.se <- SE_Filt[, SE_Filt$Timepoint != 'd0']

pcaORGsGenotype1 <- pcaSE(ORG.se, PlotTitle='First-Second Component CO', condition="Genotype", Interactive=TRUE)
pcaORGsGenotype2 <- pcaSE(ORG.se, PlotTitle='Second-Third Component CO', components = c(2,3), condition="Genotype", Interactive=TRUE)

pcaORGsTimepoint1 <- pcaSE(ORG.se, PlotTitle='First-Second Component CO', condition="Timepoint", Interactive=TRUE)
pcaORGsTimepoint2 <- pcaSE(ORG.se, PlotTitle='Second-Third Component CO', components = c(2,3), condition="Timepoint", Interactive=TRUE)

pcaORGsBatch <- pcaSE(ORG.se, PlotTitle='First-Second Component CO', condition="Batch", Interactive=TRUE)

Genotype PC1-PC2

pcaORGsGenotype1

Genotype PC2-PC3

pcaORGsGenotype2

Timepoint PC1-PC2

pcaORGsTimepoint1

Timepoint PC2-PC3

pcaORGsTimepoint2

Batch PC1-PC2

pcaORGsBatch

11. Savings

Saved objects

Intial SE without any filtering
SE_Filt after expressed gene selection and normalization
Analysis image

saveRDS(SE, paste0(SavingFolder, '', Dataset, 'SE.rds')) 

saveRDS(SE_Bio, paste0(SavingFolder, '', Dataset, 'SE_Bio.rds'))

saveRDS(SE_Filt, paste0(SavingFolder, '', Dataset, 'SE_Filt.rds')) 

SessionInfo <- sessionInfo()
Date <- date()
save.image(paste0(SavingFolder, '', Dataset, 'ExploratoryAnalysisWorkspace.RData'))

SessionInfo
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] enrichplot_1.14.2           clusterProfiler_4.2.2       circlize_0.4.15            
##  [4] SEtools_1.8.0               ComplexHeatmap_2.10.0       pheatmap_1.0.12            
##  [7] RColorBrewer_1.1-3          apeglm_1.16.0               DEFormats_1.22.0           
## [10] UpSetR_1.4.0                affy_1.72.0                 ashr_2.2-54                
## [13] sechm_1.2.0                 DESeq2_1.34.0               GenomicFeatures_1.46.5     
## [16] rtracklayer_1.54.0          tidyr_1.2.0                 biomaRt_2.50.3             
## [19] RNASeqBulkExploratory_0.2.1 org.Hs.eg.db_3.14.0         AnnotationDbi_1.56.2       
## [22] edgeR_3.36.0                limma_3.50.1                dplyr_1.0.9                
## [25] ggplot2_3.5.1               SummarizedExperiment_1.24.0 Biobase_2.54.0             
## [28] GenomicRanges_1.46.1        GenomeInfoDb_1.30.1         IRanges_2.28.0             
## [31] S4Vectors_0.32.4            BiocGenerics_0.40.0         MatrixGenerics_1.6.0       
## [34] matrixStats_0.62.0          gridExtra_2.3               DT_0.25                    
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3           SparseM_1.81             coda_0.19-4             
##   [4] bit64_4.0.5              knitr_1.40               ggvenn_0.1.10           
##   [7] irlba_2.3.5.1            DelayedArray_0.20.0      data.table_1.14.2       
##  [10] KEGGREST_1.34.0          RCurl_1.98-1.7           doParallel_1.0.17       
##  [13] generics_0.1.3           preprocessCore_1.56.0    RSQLite_2.2.12          
##  [16] shadowtext_0.1.2         bit_4.0.4                xml2_1.3.3              
##  [19] assertthat_0.2.1         viridis_0.6.2            xfun_0.31               
##  [22] hms_1.1.2                jquerylib_0.1.4          evaluate_0.17           
##  [25] TSP_1.2-1                fansi_0.5.0              restfulr_0.0.15         
##  [28] progress_1.2.2           dbplyr_2.2.0             igraph_1.3.5            
##  [31] DBI_1.1.3                geneplotter_1.72.0       htmlwidgets_1.5.4       
##  [34] purrr_0.3.5              ellipsis_0.3.2           crosstalk_1.2.0         
##  [37] backports_1.4.1          V8_4.2.0                 annotate_1.72.0         
##  [40] vctrs_0.6.5              cachem_1.0.6             withr_2.5.0             
##  [43] ggforce_0.4.1            BSgenome_1.62.0          bdsmatrix_1.3-6         
##  [46] checkmate_2.1.0          treeio_1.18.1            GenomicAlignments_1.30.0
##  [49] prettyunits_1.1.1        cluster_2.1.2            DOSE_3.20.1             
##  [52] ape_5.6-2                lazyeval_0.2.2           crayon_1.4.1            
##  [55] genefilter_1.76.0        pkgconfig_2.0.3          labeling_0.4.2          
##  [58] tweenr_2.0.2             nlme_3.1-152             seriation_1.3.6         
##  [61] rlang_1.1.4              lifecycle_1.0.3          downloader_0.4          
##  [64] registry_0.5-1           filelock_1.0.2           affyio_1.64.0           
##  [67] BiocFileCache_2.2.1      invgamma_1.1             polyclip_1.10-0         
##  [70] graph_1.72.0             Matrix_1.3-4             aplot_0.1.8             
##  [73] GlobalOptions_0.1.2      png_0.1-7                viridisLite_0.4.1       
##  [76] rjson_0.2.21             bitops_1.0-7             Biostrings_2.62.0       
##  [79] blob_1.2.3               shape_1.4.6              mixsqp_0.3-43           
##  [82] stringr_1.4.1            SQUAREM_2021.1           qvalue_2.26.0           
##  [85] gridGraphics_0.5-1       scales_1.3.0             memoise_2.0.1           
##  [88] magrittr_2.0.3           plyr_1.8.7               zlibbioc_1.40.0         
##  [91] scatterpie_0.1.7         compiler_4.1.1           tinytex_0.42            
##  [94] BiocIO_1.4.0             bbmle_1.0.25             clue_0.3-61             
##  [97] Rsamtools_2.10.0         cli_3.6.3                XVector_0.34.0          
## [100] patchwork_1.1.2          MASS_7.3-54              mgcv_1.8-37             
## [103] tidyselect_1.2.0         stringi_1.7.8            highr_0.9               
## [106] emdbook_1.3.12           yaml_2.3.5               GOSemSim_2.20.0         
## [109] locfit_1.5-9.6           ggrepel_0.9.1            sass_0.4.2              
## [112] fastmatch_1.1-3          randomcoloR_1.1.0.1      tools_4.1.1             
## [115] parallel_4.1.1           rstudioapi_0.14          foreach_1.5.2           
## [118] farver_2.1.1             Rtsne_0.16               ggraph_2.0.5            
## [121] digest_0.6.28            BiocManager_1.30.18      Rcpp_1.0.9              
## [124] httr_1.4.4               colorspace_2.0-3         XML_3.99-0.10           
## [127] topGO_2.46.0             truncnorm_1.0-8          splines_4.1.1           
## [130] yulab.utils_0.0.5        tidytree_0.3.9           graphlayouts_0.8.2      
## [133] ggplotify_0.1.0          plotly_4.10.4            xtable_1.8-4            
## [136] jsonlite_1.7.2           ggtree_3.2.1             tidygraph_1.2.1         
## [139] ggfun_0.0.7              R6_2.5.1                 pillar_1.8.1            
## [142] htmltools_0.5.2          glue_1.6.2               fastmap_1.1.0           
## [145] BiocParallel_1.28.3      codetools_0.2-18         fgsea_1.20.0            
## [148] mvtnorm_1.1-3            utf8_1.2.2               lattice_0.20-45         
## [151] bslib_0.4.0              tibble_3.1.7             sva_3.42.0              
## [154] numDeriv_2016.8-1.1      curl_4.3.3               zip_2.2.1               
## [157] GO.db_3.14.0             openxlsx_4.2.5           survival_3.2-13         
## [160] rmarkdown_2.17           munsell_0.5.0            DO.db_2.9               
## [163] GetoptLong_1.0.5         GenomeInfoDbData_1.2.7   iterators_1.0.14        
## [166] reshape2_1.4.4           gtable_0.3.1
Date
## [1] "Fri Jan 10 18:42:20 2025"

Bulk RNASeq of CHD2 iPSCs and Organoids: Exploratory Analysis

Oliviero Leonardi

January 10, 2025

All samples: iPSCs, Day25 CO, Day50 CO

2. Data Upload

2.1 Read Count Matrix

2.2 Retrieve Annotation from GeneAnnotation table and match Gene columns

2.3 GeneAnnotation table sanity checks

2.4 Specify experimental design and generate sample table

2.4.1 Import SampleSheet

2.4.2 Interactive table

2.5 Set ENSG as rownames and drop gene columns (but keep for later use)

3. Read Counts density

4. SummarizedExperiment

4.1 Consistency checks

4.2 Create SummatizedExperiment object (direct function)

4.3 Set ENSG with version as rownames

5. Biotypes selection

5.1 Biotype-selected SE

5.2 Diagnostic Barplots

Percentage of Reads associated to most expressed gene

Percentage of Removed Reads after biotype selection

6. Gene expression filtering and SE normalization

6.1 Setting gene filtering thresholds

6.2 Filter out low expressed genes

6.3 Normalize

7. Diagnostic barplots on normalized data

7.1 Barplot for library size

Library size before filtering

Library size after filtering

7.2 Barplot for TMM and expressed genes

TMM Barplot

Barplot for genes with Cpm higher than threshold

8. Density plot for Count distribution after filtering

9. Sample-to-sample correlation

10. Principal Component Analysis

10.1 PCA with Features

All Samples

Organoids

Organoids day25

Organoids day50

10.2 PCA minimal

Genotype PC1-PC2

Genotype PC2-PC3

Timepoint PC1-PC2

Timepoint PC2-PC3

Batch PC1-PC2

10.3 Timepoints-wise

iPSCs

Day25

Day50

10.3 Organoids

Genotype PC1-PC2

Genotype PC2-PC3

Timepoint PC1-PC2

Timepoint PC2-PC3

Batch PC1-PC2

11. Savings

last update on: Fri Jan 10 18:44:22 2025